Paqueterias

Definimos los paquetes que se utilizarán

## Repositorio CRAN
cran_packages <- c("bookdown", "knitr", "tidyverse", "plyr", "grid",
                   "gridExtra", "kableExtra", "xtable", "ggpubr")
## Repositorio Bioconductor
bioc_packages <- c("phyloseq", "dada2", "DECIPHER", "phangorn", "ggpubr",
                   "BiocManager","DESeq2", "microbiome", "philr")
## Repositorio GitHub
git_source <- c("twbattaglia/btools", "gmteunisse/fantaxtic",
                "MadsAlbertsen/ampvis2", "opisthokonta/tsnemicrobiota")
# fuente/nombre del paquete
git_packages <- c("btools", "fantaxtic", "ampvis2", "tsnemicrobiota")

Instalamos los paquetes, para ello es necesario que descomentes las siguientes líneas.

# # Intalar paquetes CRAN
# .inst <- cran_packages %in% installed.packages()
# if(any(!.inst)) {
#   install.packages(cran_packages[!.inst])
# }
# # Intalar paquetes BioConductor
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# .inst <- bioc_packages %in% installed.packages()
# if(any(!.inst)) {
#   BiocManager::install(bioc_packages[!.inst])
# }
# # Instalar paquetes GitHub
# .inst <- git_source %in% installed.packages()
# install.packages("devtools")
# if(any(!.inst)) {
#   devtools::install_github(git_source[!.inst])
# }
# library("devtools")
# install_github("opisthokonta/tsnemicrobiota")
# install.packages("remotes")
# remotes::install_github("kasperskytte/ampvis2")
# devtools::install_github('twbattaglia/btools')

Cargamos los paquetes

sapply(c(cran_packages, bioc_packages, git_packages), require, 
       character.only = TRUE)

##       bookdown          knitr      tidyverse           plyr           grid 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##      gridExtra     kableExtra         xtable         ggpubr       phyloseq 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##          dada2       DECIPHER       phangorn         ggpubr    BiocManager 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##         DESeq2     microbiome          philr         btools      fantaxtic 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##        ampvis2 tsnemicrobiota 
##           TRUE           TRUE

library("phyloseq")
library("ggplot2")
library("readr")
library("patchwork")
library("pheatmap")
library("microbiome")
library("xtable")
library("tsnemicrobiota")
library("Rtsne")
library("ampvis2")
library("btools")
library("RColorBrewer")
library("igraph")
library("vegan")
library("factoextra")
library("pbkrtest")
library("BiodiversityR")

Directorio de trabajo

Ajustamos directorio de trabajo donde se encuentran las lecturas “reads”

# IMPORTANTE
# Ajustamos path de trabajo según tu PC
# Mostramos directorio actual
getwd()

## [1] "C:/Users/alber/OneDrive - CINVESTAV/Documentos/Doctorado/ProyectosR/Calakmul/Shotgun"

# Si es necesario, cambiar la ruta donde están los archivos almacenados. 
# "." significa el directorio actual. 
workingDir <-  "."
setwd(workingDir)
# Creamos directorio de Resultados y otras capertas
dir.create("Resultados/")

## Warning in dir.create("Resultados/"): 'Resultados' already exists

dir.create("Resultados/Abundancias/")

## Warning in dir.create("Resultados/Abundancias/"): 'Resultados\Abundancias'
## already exists

dir.create("Resultados/AbundanciasAV2/")

## Warning in dir.create("Resultados/AbundanciasAV2/"): 'Resultados\AbundanciasAV2'
## already exists

dir.create("Resultados/BetaDiversidad/")

## Warning in dir.create("Resultados/BetaDiversidad/"): 'Resultados\BetaDiversidad'
## already exists

# Asignamos la carpeta donde están las lecturas crudas
miseq_path <- "RawData"
list.files(miseq_path)

## character(0)

Importación y manipulación de archivo kraken-biom

Metemos a la variable merged_metagenomes el resultado de kraken-biom, esto resulta en un objeto phyloseq

merged_metagenomes <- import_biom("CalakmulShotGun.biom")
class(merged_metagenomes)

## [1] "phyloseq"
## attr(,"package")
## [1] "phyloseq"

head(merged_metagenomes@tax_table@.Data)

##         Rank1         Rank2               Rank3                   
## 356     "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 41294   "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 374     "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 1437360 "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 1274631 "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 2057741 "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
##         Rank4            Rank5                  Rank6              
## 356     "o__Rhizobiales" "f__"                  "g__"              
## 41294   "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__"              
## 374     "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
## 1437360 "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
## 1274631 "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
## 2057741 "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
##         Rank7            
## 356     "s__"            
## 41294   "s__"            
## 374     "s__"            
## 1437360 "s__erythrophlei"
## 1274631 "s__icense"      
## 2057741 "s__sp. SK17"

En los datos, el nombre como tal empieza en la letra numero 4, tiene 3 antes, por ejemplo “k__“, con el siguiente comando eliminamos estas 3 primeras letras y empieza en la 4ta.

merged_metagenomes@tax_table@.Data <- substring(merged_metagenomes@tax_table@.Data, 4)
# Cambiamos el nombre de las columnas con el siguiente comando
colnames(merged_metagenomes@tax_table@.Data) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
head(merged_metagenomes@tax_table@.Data)

##         Kingdom    Phylum           Class                 Order        
## 356     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 41294   "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 374     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1437360 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1274631 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 2057741 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
##         Family              Genus            Species       
## 356     ""                  ""               ""            
## 41294   "Bradyrhizobiaceae" ""               ""            
## 374     "Bradyrhizobiaceae" "Bradyrhizobium" ""            
## 1437360 "Bradyrhizobiaceae" "Bradyrhizobium" "erythrophlei"
## 1274631 "Bradyrhizobiaceae" "Bradyrhizobium" "icense"      
## 2057741 "Bradyrhizobiaceae" "Bradyrhizobium" "sp. SK17"

# Lo metemos en una tabla
TablaMergedMetagenomes <- merged_metagenomes@tax_table@.Data
# Guardamos
write.table(TablaMergedMetagenomes, file = "CalakmulTablaMergedMetagenomes.txt", 
            row.names=TRUE, col.names=NA, quote=FALSE, sep="\t")

Podemos ver los únicos Kingdom, Phylum, Class, etc…

head(unique(merged_metagenomes@tax_table@.Data[,"Kingdom"]))

## [1] "Bacteria"

head(unique(merged_metagenomes@tax_table@.Data[,"Phylum"]))

## [1] "Proteobacteria"      "Actinobacteria"      "Firmicutes"         
## [4] "Deinococcus-Thermus" "Cyanobacteria"       "Chloroflexi"

head(unique(merged_metagenomes@tax_table@.Data[,"Class"]))

## [1] "Alphaproteobacteria"   "Betaproteobacteria"    "Gammaproteobacteria"  
## [4] "Deltaproteobacteria"   "Epsilonproteobacteria" "Acidithiobacillia"

head(unique(merged_metagenomes@tax_table@.Data[,"Order"]))

## [1] "Rhizobiales"      "Sphingomonadales" "Rhodospirillales" "Rhodobacterales" 
## [5] "Caulobacterales"  "Rickettsiales"

head(unique(merged_metagenomes@tax_table@.Data[,"Family"]))

## [1] ""                    "Bradyrhizobiaceae"   "Rhizobiaceae"       
## [4] "Methylobacteriaceae" "Hyphomicrobiaceae"   "Phyllobacteriaceae"

head(unique(merged_metagenomes@tax_table@.Data[,"Genus"]))

## [1] ""                 "Bradyrhizobium"   "Rhodopseudomonas" "Bosea"           
## [5] "Afipia"           "Nitrobacter"

head(unique(merged_metagenomes@tax_table@.Data[,"Species"]))

## [1] ""               "erythrophlei"   "icense"         "sp. SK17"      
## [5] "lablabi"        "diazoefficiens"

Como en Kingdom hay “solo” bacterias, nos vamos a quedar solo con eso, y así quitar lo que no entra en esta categoria.

merged_metagenomes <- subset_taxa(merged_metagenomes, Kingdom == "Bacteria")

Para ver cuantos de un cierto taxón hay dentro de otro

# Aquí vemos cuantos Firmicutes hay en el Phylum
sum(merged_metagenomes@tax_table@.Data[,"Phylum"] == "Firmicutes")

## [1] 708

Añadimos metadatos

MetaShotgun <- data.frame(SAMPLE=c("Ag1", "Ag2", "Ag3"), 
                          SITE=c("2", "3", "1"), 
                          pH=c(6.2, 6.7, 6.1), 
                          N_TOTAL=c(0.36, 0.22, 0.29),
                          P_OLSEN=c(29.6, 13.5, 4.0),
                          ORGANIC=c(12.0, 13.5, 24.2),
                          SITE_ID=c("Ag1", "Ag2", "Ag3"))
rownames(MetaShotgun) <- MetaShotgun$SAMPLE
rownames(merged_metagenomes) %in% rownames(MetaShotgun)

## logical(0)

all(rownames(merged_metagenomes) %in% MetaShotgun$SAMPLE)

## [1] TRUE

merged_metagenomes@sam_data <- sample_data(MetaShotgun)

Ensamble y manipulación de objetos de Phyloseq

Resumen del objeto Phyloseq

# Nos dice el número de taxas en las muestras y en cuantos rangos taxonómicos
merged_metagenomes

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 4457 taxa and 3 samples ]
## sample_data() Sample Data:       [ 3 samples by 7 sample variables ]
## tax_table()   Taxonomy Table:    [ 4457 taxa by 7 taxonomic ranks ]

# Cantidad de elementos en cada muestra
sample_sums(merged_metagenomes)

##     Ag1     Ag2     Ag3 
## 1633465 1844772 1611481

# Un resumen de cada muestra, indicando cuantos mínimos medios y máximos hay para cada taxa
summary(merged_metagenomes@otu_table@.Data)

##       Ag1               Ag2               Ag3         
##  Min.   :    0.0   Min.   :    0.0   Min.   :    0.0  
##  1st Qu.:    8.0   1st Qu.:    9.0   1st Qu.:    8.0  
##  Median :   39.0   Median :   38.0   Median :   37.0  
##  Mean   :  366.5   Mean   :  413.9   Mean   :  361.6  
##  3rd Qu.:  307.0   3rd Qu.:  328.0   3rd Qu.:  295.0  
##  Max.   :40397.0   Max.   :49917.0   Max.   :39700.0

Gráfico

# Generamos un data.frame con datos del merge_metagenomes
deep <- data.frame(Samples = sample_names(merged_metagenomes),
                    Reads = sample_sums(merged_metagenomes)/1000000)
deep

##     Samples    Reads
## Ag1     Ag1 1.633465
## Ag2     Ag2 1.844772
## Ag3     Ag3 1.611481

# Graficamos
COLOR = c("#a6cee3", "#b2df8a", "#fdbf6f")
plot.reads <- ggplot(data = deep, mapping = aes(x = Samples, y = Reads, fill = Samples)) +
  theme_classic() +
  labs(x="Samples", y="Million reads",
       title = "Total readings of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_col(position = position_dodge(0.9), width = 0.9) +
  scale_fill_brewer(palette = "Set1")
# Verla en el markdown
plot.reads

pdf("PlotReads.pdf", width=10, height=7)
plot.reads
dev.off()

## png 
##   2

svg("PlotReads.svg", width=10, height=7)
plot.reads
dev.off()

## png 
##   2

Estructura y manipulación de un objeto phyloseq

Muchas veces queremos analizar un sub conjunto de las muestras en nuestro objeto phyloseq, o bien, queremos seleccionar ciertos grupos taxonómicos para análisis posteriores. Phyloseq nos permite hacer todo tipo de filtros para llevar esto a cabo. Veamos dónde se almacena la información en phyloseq.

# Número de taxas
ntaxa(merged_metagenomes)

## [1] 4457

# Número de muestras
nsamples(merged_metagenomes)

## [1] 3

# Nombre de las muestras
sample_names(merged_metagenomes)

## [1] "Ag1" "Ag2" "Ag3"

# Niveles taxonómicos 
rank_names(merged_metagenomes)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

# Nombre de las variables de los metadatos
sample_variables(merged_metagenomes)

## [1] "SAMPLE"  "SITE"    "pH"      "N_TOTAL" "P_OLSEN" "ORGANIC" "SITE_ID"

También podemos ver un resumen estadístico de la tabla de OTUs

# Resumen estadístico
summary(merged_metagenomes@otu_table@.Data)

##       Ag1               Ag2               Ag3         
##  Min.   :    0.0   Min.   :    0.0   Min.   :    0.0  
##  1st Qu.:    8.0   1st Qu.:    9.0   1st Qu.:    8.0  
##  Median :   39.0   Median :   38.0   Median :   37.0  
##  Mean   :  366.5   Mean   :  413.9   Mean   :  361.6  
##  3rd Qu.:  307.0   3rd Qu.:  328.0   3rd Qu.:  295.0  
##  Max.   :40397.0   Max.   :49917.0   Max.   :39700.0

La tabla de taxonomía relaciona el nombre de las taxa con el linaje taxonómico de éstas, por ejemplo, vincula una variante de secuencia, o ASV, con los rangos taxonómicos desde Reino hasta Género o Especie dependiendo del nivel de resolución del análisis.

# Tabla de taxonomía
tax_table(merged_metagenomes)[1:5, 1:4]

## Taxonomy Table:     [5 taxa by 4 taxonomic ranks]:
##         Kingdom    Phylum           Class                 Order        
## 356     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 41294   "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 374     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1437360 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1274631 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"

Ahora, el objeto phyloseq se ha vuelto un estándar en la industria ya que otros paquetes ahora usan esta estructura de datos para sus propias funciones. Uno de esos paquetes es microbiome y ampvis. Podemos fácilmente obtener un resumen global de nuestro objeto phyloseq usando la función summarize_phyloseq.

summarize_phyloseq(merged_metagenomes)

## [[1]]
## [1] "1] Min. number of reads = 1611481"
## 
## [[2]]
## [1] "2] Max. number of reads = 1844772"
## 
## [[3]]
## [1] "3] Total number of reads = 5089718"
## 
## [[4]]
## [1] "4] Average number of reads = 1696572.66666667"
## 
## [[5]]
## [1] "5] Median number of reads = 1633465"
## 
## [[6]]
## [1] "7] Sparsity = 0.0270735173135891"
## 
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
## 
## [[8]]
## [1] "8] Number of singletons = 65"
## 
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n        (i.e. exactly one read detected across all samples)1.4583800762845"
## 
## [[10]]
## [1] "10] Number of sample variables are: 7"
## 
## [[11]]
## [1] "SAMPLE"  "SITE"    "pH"      "N_TOTAL" "P_OLSEN" "ORGANIC" "SITE_ID"

Este comando nos muestra el mínimo y máximo de reads, número total y promedio de reads, etc. También muestra los encabezados de las columans en la tabla de metadata. Veamos ahora una tabla que mezcla metadata, taxonomía y abundancia del taxon más abundante de cada muestra.

df <- psmelt(merged_metagenomes)
head(df)

##          OTU Sample Abundance SAMPLE SITE  pH N_TOTAL P_OLSEN ORGANIC SITE_ID
## 4800    1883    Ag2     49917    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## 4798    1883    Ag1     40397    Ag1    2 6.2    0.36    29.6    12.0     Ag1
## 4799    1883    Ag3     39700    Ag3    1 6.1    0.29     4.0    24.2     Ag3
## 5126  191495    Ag2     36610    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## 9155     374    Ag2     25994    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## 11663 674703    Ag2     25729    Ag2    3 6.7    0.22    13.5    13.5     Ag2
##        Kingdom         Phylum               Class               Order
## 4800  Bacteria Actinobacteria      Actinobacteria    Streptomycetales
## 4798  Bacteria Actinobacteria      Actinobacteria    Streptomycetales
## 4799  Bacteria Actinobacteria      Actinobacteria    Streptomycetales
## 5126  Bacteria Actinobacteria     Thermoleophilia Solirubrobacterales
## 9155  Bacteria Proteobacteria Alphaproteobacteria         Rhizobiales
## 11663 Bacteria Proteobacteria Alphaproteobacteria         Rhizobiales
##                  Family          Genus       Species
## 4800  Streptomycetaceae   Streptomyces              
## 4798  Streptomycetaceae   Streptomyces              
## 4799  Streptomycetaceae   Streptomyces              
## 5126  Conexibacteraceae   Conexibacter        woesei
## 9155  Bradyrhizobiaceae Bradyrhizobium              
## 11663 Hyphomicrobiaceae    Rhodoplanes sp. Z2-YC6860

Ahora veamos cómo podemos filtrar y hacer subsetting de un objeto phyloseq. Esto lo hacemos con tres grupos de funciones, por ejemplo, filter, subset, y prune. + Filtrar se refiere a filtrar según alguna regla lógica. Ya lo hicimos en la parte de control de calidad cuando llamamos la función filter_taxa(psd1, function(x) mean(x) > 1e-5, TRUE). Acá le pedíamos a la función filter_taxa sobre el objeto merged_metagenomes, calculara la media de cuentas de lecturas para cada taxa y si este resultado era menor que 1e-5, lo eliminara. Ahora filtramos según abundancia. Primero transformamos en abundancia relativa y luego filtramos.

# Transformamos las cuentas en porcentaje
Porcentaje <- transform_sample_counts(merged_metagenomes, function(x) x / sum(x) * 100 )
# Filtramos las taxa con una abundancia inferior al 1%
Porcentaje.filtrado <- filter_taxa(Porcentaje, function(x) sum(x) > 1, TRUE)

¿Cuántas taxa permanecen en nuestro objeto phyloseq? Con una operación tan simple como la que acabamos de aplicar, nos damos cuenta que la mayoría de las taxa presentes en nuestras están en muy baja abundancia. Ahora imaginemos la situación donde queremos filtrar nuestro objeto pero en función de un taxon en específico.

# Ahora filtramos de acuerdo a Acidobacteriota
Porcentaje.filtrado.Proteobacteria <- subset_taxa(Porcentaje.filtrado, Phylum == "Proteobacteria")
# También podríamos todo lo que NO es Acidobacteriota
Porcentaje.filtrado.NoAcidobacteriota <- subset_taxa(Porcentaje.filtrado, Phylum != "Acidobacteriota")

Otra manera de filtrar un objeto phyloseq es en base a algún atributo presente en sample_data. Por ejemplo, con estos datos uno podría querer estudiar el microbioma de área por separado. Para esto crearíamos tres objetos phyloseq a partir de merged_metagenomes.

merged_metagenomes.Ag1 <- subset_samples(merged_metagenomes, SITE_ID == "Ag1")
merged_metagenomes.Ag2 <- subset_samples(merged_metagenomes, SITE_ID == "Ag2")
merged_metagenomes.Ag3 <- subset_samples(merged_metagenomes, SITE_ID == "Ag3")

Alternativamente, podríamos decidir estudiar solo dos de las 3 áreas

merged_metagenomes.Ag1_Ag2 <- subset_samples(subset_samples(merged_metagenomes, SITE_ID != "Ag3"))

Transformación y manipulación de datos

Siguimos viendo el objeto phyloseq

# Con esto nos dice cuántos datos vacios hay, en este caso serán todos aquellos que aparezcan como TRUE
summary(merged_metagenomes@tax_table@.Data== "") # Hay 1337 NAs

##   Kingdom          Phylum          Class           Order        
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4457      FALSE:4457      FALSE:4313      FALSE:4445     
##                                  TRUE :144       TRUE :12       
##    Family          Genus          Species       
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4381      FALSE:4217      FALSE:3807     
##  TRUE :76        TRUE :240       TRUE :650

#Actualizamos la variable merged_metagenomes indicando que se incluya todo lo diferente (!=) a vacio "", para la variable "Genus"
merged_metagenomes <- subset_taxa(merged_metagenomes, Genus != "")
summary(merged_metagenomes@tax_table@.Data== "") # Ya no hay ninguno

##   Kingdom          Phylum          Class           Order        
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4217      FALSE:4217      FALSE:4088      FALSE:4210     
##                                  TRUE :129       TRUE :7        
##    Family          Genus          Species       
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4206      FALSE:4217      FALSE:3791     
##  TRUE :11                        TRUE :426

Datos Absolutos_S y Relativos_S

Generamos una tabla de datos Absolutos_S y Relativos_S a nivel de “Phylum”, para posteriores análisis.

Relativos_S <- transform_sample_counts(merged_metagenomes, function(x) x*100 / sum(x))
# Valores Relativos_S
Glom_Rel_S <- tax_glom(physeq = Relativos_S, taxrank = 'Phylum')
# psmelt permite generar un data.frame a partir de un objeto de phyloseq
Relativos_S <- psmelt(Glom_Rel_S)
# Valores Absolutos_S
Glom_Abs_S <- tax_glom(physeq = merged_metagenomes, taxrank = "Phylum")
Absolutos_S <- psmelt(Glom_Abs_S)

Gráficas de valores Absolutos_S y Relativos_S

En eje X usamos la variable “SITE_ID”, y en Y es “Abundance” y rellenamos por “Phylum”

Abs.plot <- ggplot(data=Absolutos_S, aes(x=SITE_ID, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="stack") + 
  theme_classic() +
  labs(x="Samples", y="Absolute abundance",
       title = "Absolute abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))

Rel.plot <- ggplot(data=Relativos_S, aes(x=SITE_ID, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="fill") + 
  theme_classic() +
  labs(x="Samples", y="Relative abundance",
       title = "Relative abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))

Abs.plot | Rel.plot

Abundancia absoluta y relativa a nivel de Phylum.

pdf("Resultados/Abundancias/PlotRaw.pdf", width=10, height=7)
Abs.plot
dev.off()

## png 
##   2

svg("Resultados/Abundancias/PlotRaw.svg", width=10, height=7)
Abs.plot
dev.off()

## png 
##   2

pdf("Resultados/Abundancias/PlotRel.pdf", width=10, height=7)
Rel.plot
dev.off()

## png 
##   2

svg("Resultados/Abundancias/PlotRel.svg", width=10, height=7)
Rel.plot
dev.off()

## png 
##   2

pdf("Resultados/Abundancias/PlotRaw-Rel.pdf", width=15, height=7)
Abs.plot | Rel.plot
dev.off()

## png 
##   2

svg("Resultados/Abundancias/PlotRaw-Rel.svg", width=15, height=7)
Abs.plot | Rel.plot
dev.off()

## png 
##   2

Podemos hacerlo para las diferentes variables en X, Y y fill

names(Relativos_S)

##  [1] "OTU"       "Sample"    "Abundance" "SAMPLE"    "SITE"      "pH"       
##  [7] "N_TOTAL"   "P_OLSEN"   "ORGANIC"   "SITE_ID"   "Kingdom"   "Phylum"

names(Absolutos_S)

##  [1] "OTU"       "Sample"    "Abundance" "SAMPLE"    "SITE"      "pH"       
##  [7] "N_TOTAL"   "P_OLSEN"   "ORGANIC"   "SITE_ID"   "Kingdom"   "Phylum"

Renombramos aquellos “Phylum” que presenten una abundancia menor a 1 %.

Relativos_S$Phylum <- as.character(Relativos_S$Phylum) # Primero lo convertimos a letra
Relativos_S$Phylum[Relativos_S$Abundance < 1.0] <- "Phylum < 1.0% abund." # Renombramos
unique(Relativos_S$Phylum)

## [1] "Proteobacteria"       "Actinobacteria"       "Firmicutes"          
## [4] "Planctomycetes"       "Acidobacteria"        "Phylum < 1.0% abund."

Volvemos a graficar usando ese filtro de “Phylum < 1.0% abund.”

Rel.plot2 <- ggplot(data=Relativos_S, aes(x=SITE_ID, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="fill") + 
  theme_classic() +
  labs(x="Samples", y="Relative abundance",
       title = "Relative abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))

Abs.plot | Rel.plot2

Abundancia absoluta y relativa a nivel de Phylum. con Phylum < 1.0 % abundancia.

pdf("Resultados/Abundancias/PlotRel2.pdf", width=10, height=7)
Rel.plot2
dev.off()

## png 
##   2

svg("Resultados/Abundancias/PlotRel2.svg", width=10, height=7)
Rel.plot2
dev.off()

## png 
##   2

pdf("Resultados/Abundancias/PlotRaw-Rel2.pdf", width=15, height=7)
Abs.plot | Rel.plot2
dev.off()

## png 
##   2

svg("Resultados/Abundancias/PlotRaw-Rel2.svg", width=15, height=7)
Abs.plot | Rel.plot2
dev.off()

## png 
##   2

Abundancias

Graficamos con ggplot2 para ver el comportamiento de las lecturas en cada una de las muestras

#Generamos un data.frame con datos del merge_metagenomes
deep <- data.frame(Samples = merged_metagenomes@sam_data@.Data[[1]],
                   Reads = sample_sums(merged_metagenomes))
deep

##     Samples   Reads
## Ag1     Ag1 1523762
## Ag2     Ag2 1726131
## Ag3     Ag3 1503461

#Graficamos
plot.reads <- ggplot(data = deep, mapping = aes(x = Samples, y = Reads, fill = Samples)) +
  theme_classic() +
  labs(x="Samples", y="Reads",
       title = "Total readings of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_col(position = position_dodge(0.9), width = 0.9) +
  scale_fill_brewer(palette = "Paired")
plot.reads

Distribución del total de lecturas en cada una de las muestras.

pdf("Resultados/Abundancias/PlotReads.pdf", width=10, height=7)
plot.reads
dev.off()

## png 
##   2

svg("Resultados/Abundancias/PlotReads.svg", width=10, height=7)
plot.reads
dev.off()

## png 
##   2

Nos quedamos con el nombre de aquellos Phylum con procentaje mayor a 1

OnePercentage <- unique(Relativos_S$Phylum[Relativos_S$Abundance > 1.0])
OnePercentage

## [1] "Proteobacteria" "Actinobacteria" "Firmicutes"     "Planctomycetes"
## [5] "Acidobacteria"

Abundancias a nivel de Phylum nuevamente para los de abundancia mayor al 1 %

for (i in OnePercentage){
  tryCatch({
    # Prueba 
    # i <- "Actinobacteriota"
    phy <- subset_taxa(merged_metagenomes, Kingdom == "Bacteria") # Nos quedamos solo con Bacterias
    phy <- subset_taxa(merged_metagenomes, Phylum == i) # Nos quedamos con el phylum de interés
    # Comprobamos
    paste("El siguiente gráfico es del Phylum", unique(phy@tax_table@.Data[,2])) 
    # Lo hacemos valores Relativos_S
    phy <- transform_sample_counts(phy, function(x) x*100 / sum(x) ) 
    # Cortamos a nivel de género
    phy <- tax_glom(phy, taxrank = "Genus")
    # Lo volvemos data.frame
    phy.t <- psmelt(phy)
    # Renombramos los géneros que tienen abundancia menor a 1 %
    phy.t$Genus[phy.t$Abundance < 1.0] <- "Genus < 1.0 % abundance"
    }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) 
    # Graficamos resécto a "SITE_ID"
    plot <- ggplot(data = phy.t, aes(x=SITE_ID, y=Abundance, fill=Genus)) + 
      geom_bar(aes(), stat="identity", position="fill") +
      theme(legend.position = "right", 
            legend.text = element_text(size=12, face="bold"),
            legend.title = element_text(size = 12, face = "bold"),
            axis.text = element_text(color = "black", size = 12, face = "bold"),
            axis.title = element_text(color = "black", size = 12, face = "bold"),
            plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
      labs(x="Samples", y="Relative abundance",
          title = paste("Abundancia relativa dentro del phylum", i))
    print(plot)
    pdf(paste0("Resultados/Abundancias/PlotRel", i, ".pdf"), width = 30, height = 10)
    print(plot)
    dev.off()
    svg(paste0("Resultados/Abundancias/PlotRel", i, ".svg"), width = 30, height = 10)
    print(plot)
    dev.off()
}

Abundancia relativa de todos los Phylum presentes en mayor a 1 % de Abundancia.

Análisis de diversidad

¿Qué entendemos por diversidad? Al menos podemos conceptualizar diversidad a dos niveles: diversidad genética o morfológica, y biodiversidad. En el estudio de comunidades, tomamos prestado el concepto de biodiversidad de ecología de comunidades donde estamos interesados en la riqueza de especies (número de especies diferentes en una comunidad o diversidad alfa), en las diferencias y similitudes entre comunidades (diversidad beta), y en algunos casos en la diversidad total de un región o paisaje ecológico (landscape; diversidad gama). Medidas de riqueza, uniformidad, dominancia, diversidad filogenética (diversidad alfa)

Alfa diversidad

En el contexto metagenómico, medimos diversidad alfa \(\alpha\) usando una serie de medidas prestadas de ecología que nos permiten caracterizar una comunidad microbiana. phyloseq tiene una función muy útil que nos permite calcular y graficar hasta siete medidas, i .e., Observed (simplemente el número de taxa o riqueza), Chao1 (la riqueza ajustada por probabilidad de no observar especies), ACE (riqueza que toma en cuenta la abundancia relativa), Shannon (abundancia relativa de taxa), Simpson (1 - la probabilidad de que observemos aleatoriamente dos bacterias en una comunidad y que pertenezcan a diferentes especies ), Inverse Simpson ( 1 / Simpson), y Fisher (riqueza tomando en cuenta abundancia). En phyloseq simplemente llamamos la función plot_richness y podemos visualizar las medidas de diversidad.

alfa.diver <- plot_richness(physeq = merged_metagenomes, 
              color =  "SITE_ID",
              x = "SITE_ID",
              measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "Fisher"),
              title = "Alpha diversity indices for the Calakmul samples obtained by ShotGun") + 
  labs(x="Samples", y="Alpha Diversity Measure") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_boxplot(aes(fill = SITE_ID), alpha=.7) + 
  scale_color_brewer(palette = "Set1") + 
  scale_fill_brewer(palette = "Set1")
alfa.diver

Índices de diversidad alpha \(lpha\).

pdf("Resultados/PlotAlfa.pdf", width=20, height=10)
alfa.diver
dev.off()

## png 
##   2

svg("Resultados/PlotAlfa.svg", width=20, height=10)
alfa.diver
dev.off()

## png 
##   2

El paquete microbiome ofrece otras herramientas para evaluar diversidad que son accesibles fácilmente a través de su función global.

Alpha <- alpha(merged_metagenomes, index = "all")
head(Alpha)

##     observed    chao1 diversity_inverse_simpson diversity_gini_simpson
## Ag1     4122 4214.560                  352.4709              0.9971629
## Ag2     4131 4214.379                  306.7775              0.9967403
## Ag3     4100 4189.446                  345.0244              0.9971017
##     diversity_shannon diversity_fisher diversity_coverage evenness_camargo
## Ag1          6.885202         515.8116                231        0.2035280
## Ag2          6.811040         508.0487                212        0.1921895
## Ag3          6.869876         513.6507                228        0.2018350
##     evenness_pielou evenness_simpson evenness_evar evenness_bulla dominance_dbp
## Ag1       0.8271414       0.08550968     0.1350859      0.3609202    0.02651136
## Ag2       0.8180177       0.07426228     0.1312996      0.3506728    0.02891843
## Ag3       0.8258311       0.08415229     0.1360249      0.3592546    0.02640574
##     dominance_dmn dominance_absolute dominance_relative dominance_simpson
## Ag1    0.04183593              40397         0.02651136       0.002837114
## Ag2    0.05012771              49917         0.02891843       0.003259692
## Ag3    0.04186274              39700         0.02640574       0.002898346
##     dominance_core_abundance dominance_gini rarity_log_modulo_skewness
## Ag1                0.2763476      0.8010570                   2.054063
## Ag2                0.2906512      0.8116390                   2.052640
## Ag3                0.2796986      0.8035928                   2.053428
##     rarity_low_abundance rarity_rare_abundance
## Ag1            0.7161781             0.7006737
## Ag2            0.7008865             0.6776919
## Ag3            0.7180931             0.6992699

La función alpha nos da 22 medidas de diversidad que nos ayudan a entender la estructura de las comunidades microbianas. En general, estas medidas se dividen en riqueza, diversidad, dominancia, rareza, cobertura y uniformidad. El paquete microbiome ofrece funciones para calcular cada uno de estos aspectos de las comunidades microbianas.

# Riqueza
Riqueza <- richness(merged_metagenomes)
Riqueza

##     observed    chao1
## Ag1     4122 4214.560
## Ag2     4131 4214.379
## Ag3     4100 4189.446

# Dominancia
Dominancia <- dominance(merged_metagenomes, index = "all")
Dominancia

##            dbp        dmn absolute   relative     simpson core_abundance
## Ag1 0.02651136 0.04183593    40397 0.02651136 0.002837114      0.2763476
## Ag2 0.02891843 0.05012771    49917 0.02891843 0.003259692      0.2906512
## Ag3 0.02640574 0.04186274    39700 0.02640574 0.002898346      0.2796986
##          gini
## Ag1 0.8010570
## Ag2 0.8116390
## Ag3 0.8035928

# Rareza
Rareza <- rarity(merged_metagenomes, index = "all")
Rareza

##     log_modulo_skewness low_abundance rare_abundance
## Ag1            2.054063     0.7161781      0.7006737
## Ag2            2.052640     0.7008865      0.6776919
## Ag3            2.053428     0.7180931      0.6992699

# Cobertura
Cobertura <- coverage(merged_metagenomes, threshold = 0.5)
Cobertura

## Ag1 Ag2 Ag3 
## 231 212 228

# Desigualdad
Desigualdad <- inequality(merged_metagenomes)
Desigualdad

##       Ag1       Ag2       Ag3 
## 0.8010570 0.8116390 0.8035928

# Uniformidad
Uniformidad <- evenness(merged_metagenomes, "all")
Uniformidad

##       camargo    pielou    simpson      evar     bulla
## Ag1 0.2035280 0.8271414 0.08550968 0.1350859 0.3609202
## Ag2 0.1921895 0.8180177 0.07426228 0.1312996 0.3506728
## Ag3 0.2018350 0.8258311 0.08415229 0.1360249 0.3592546

Graficamos los resultados y calculamos la significancia estadística. Para esto usamos el paquete ggpubr que genera “publication-ready plots”, algo que siempre es deseable (ejecuta library(ggpubr)).

library("ggpubr")
# Generamos un objeto `phyloseq` sin taxa´s que tengan 0 reads
merged_metagenomes <- prune_taxa(taxa_sums(merged_metagenomes) > 0, merged_metagenomes)
# Calculamos los índices de diversidad
Alpha2 <- alpha(merged_metagenomes, index = "all")
# Y finalmente visualizamos la tabla de resultados
head(Alpha2)

##     observed    chao1 diversity_inverse_simpson diversity_gini_simpson
## Ag1     4122 4214.560                  352.4709              0.9971629
## Ag2     4131 4214.379                  306.7775              0.9967403
## Ag3     4100 4189.446                  345.0244              0.9971017
##     diversity_shannon diversity_fisher diversity_coverage evenness_camargo
## Ag1          6.885202         515.8116                231        0.2035280
## Ag2          6.811040         508.0487                212        0.1921895
## Ag3          6.869876         513.6507                228        0.2018350
##     evenness_pielou evenness_simpson evenness_evar evenness_bulla dominance_dbp
## Ag1       0.8271414       0.08550968     0.1350859      0.3609202    0.02651136
## Ag2       0.8180177       0.07426228     0.1312996      0.3506728    0.02891843
## Ag3       0.8258311       0.08415229     0.1360249      0.3592546    0.02640574
##     dominance_dmn dominance_absolute dominance_relative dominance_simpson
## Ag1    0.04183593              40397         0.02651136       0.002837114
## Ag2    0.05012771              49917         0.02891843       0.003259692
## Ag3    0.04186274              39700         0.02640574       0.002898346
##     dominance_core_abundance dominance_gini rarity_log_modulo_skewness
## Ag1                0.2763476      0.8010570                   2.054063
## Ag2                0.2906512      0.8116390                   2.052640
## Ag3                0.2796986      0.8035928                   2.053428
##     rarity_low_abundance rarity_rare_abundance
## Ag1            0.7161781             0.7006737
## Ag2            0.7008865             0.6776919
## Ag3            0.7180931             0.6992699

Metadatos con diversidad

Extraemos el metadata del objeto phyloseq.

merged_metagenomes.meta <- meta(merged_metagenomes)
head(merged_metagenomes.meta)

##     SAMPLE SITE  pH N_TOTAL P_OLSEN ORGANIC SITE_ID
## Ag1    Ag1    2 6.2    0.36    29.6    12.0     Ag1
## Ag2    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## Ag3    Ag3    1 6.1    0.29     4.0    24.2     Ag3

Le agregamos la tabla de diversidad a la metadata

names(Alpha2)

##  [1] "observed"                   "chao1"                     
##  [3] "diversity_inverse_simpson"  "diversity_gini_simpson"    
##  [5] "diversity_shannon"          "diversity_fisher"          
##  [7] "diversity_coverage"         "evenness_camargo"          
##  [9] "evenness_pielou"            "evenness_simpson"          
## [11] "evenness_evar"              "evenness_bulla"            
## [13] "dominance_dbp"              "dominance_dmn"             
## [15] "dominance_absolute"         "dominance_relative"        
## [17] "dominance_simpson"          "dominance_core_abundance"  
## [19] "dominance_gini"             "rarity_log_modulo_skewness"
## [21] "rarity_low_abundance"       "rarity_rare_abundance"

merged_metagenomes.meta$Shannon <- Alpha2$diversity_shannon 
merged_metagenomes.meta$InverseSimpson <- Alpha2$diversity_inverse_simpson
merged_metagenomes.meta$Chao1 <- Alpha2$chao1
merged_metagenomes.meta$Observed <- Alpha2$observed
merged_metagenomes.meta$Fisher <- Alpha2$diversity_fisher

En este ejercicio nos interesa comparar la diversidad entre áreas. Recordemos que tenemos datos para tres áreas: Ag1, Ag2, Ag3. Necesitamos crear una lista de comparasiones de a pares para poder visualizar y calcular significancia estadística de manera simultánea.

# Obtenemos las variables desde nuestro objeto `phyloseq`
Area <- levels(as.factor(merged_metagenomes.meta$SITE_ID))
# Creamos una lista de lo que queremos comparar
ParesArea <- combn(seq_along(Area), 2, simplify = FALSE, FUN = function(i) Area[i])
ParesArea

## [[1]]
## [1] "Ag1" "Ag2"
## 
## [[2]]
## [1] "Ag1" "Ag3"
## 
## [[3]]
## [1] "Ag2" "Ag3"

Con la función ggviolin podemos generar un gráfico de violín en un solo paso de la siguiente forma.

# Cargamos una variable de color
# COLOR <- c("#a6cee3", "#b2df8a", "#fdbf6f")
Indices <- c("Observed","Shannon","Fisher","Chao1","InverseSimpson")
# Ciclo
for (i in Indices){
  # Prueba 
  # i <- "Observed"
  p <- ggviolin(merged_metagenomes.meta, x = "SITE_ID", y = i, 
                add = "boxplot", fill = "SITE_ID") + #, palette = COLOR)
    scale_fill_brewer(palette = "Set1") +
    # Evaluamos la significancia estadística
    stat_compare_means(comparisons = ParesArea)
  print(p)
}

## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

## Warning in max(data$density): no non-missing arguments to max; returning -Inf

## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

## Warning in max(data$density): no non-missing arguments to max; returning -Inf

## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

## Warning in max(data$density): no non-missing arguments to max; returning -Inf

## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

## Warning in max(data$density): no non-missing arguments to max; returning -Inf

## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.

## Warning in max(data$density): no non-missing arguments to max; returning -Inf

## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

Beta Diversidad

La diversidad beta \(\beta\) es la relación entre la diversidad de especies regional y local. La diversidad beta como una medida del recambio de especies enfatiza demasiado el papel de las especies raras, ya que la diferencia en la composición de especies entre dos sitios o comunidades probablemente refleje la presencia o ausencia de algunas especies raras en los ensamblajes. Existen diferentes medidas de similitud (o disimilitud, i.e., 1 - similitud) disponibles que nos permiten entender las relaciones entre nuestras muestras. En general todas producen matrices de distancia comparables. El paquete phyloseq ofrece un gran número de medidas de distancia. Las más populares son UniFrac y Weighted UniFrac (medidas que consideran filogenia) y otras independientes de filogenia como: Jaccard, Manhattan, Euclidian, Bray-Curtis, Canberra, etc. Por otra parte, la matriz de distancia resultante no se usa en aislación sino que en conjunto con algún método de ordinación o escalamiento multidimensional (ordination). De nuevo, phyloseq ofrece un gran número de métodos entre los cuales se encuentran: detrended y canonical correspondence analysis, Double Principal Coordinate Analysis, Non-metric MultiDimenstional Scaling, y MDS/PCoA. Probemos entonces hacer un análisis tipo PCoA con una matriz de distancia que considera las relaciones filogenéticas y otra que no.

# Con esto vemos todos los métodos (métricas de distancia) para calcular la beta diversidad
distanceMethodList

## $UniFrac
## [1] "unifrac"  "wunifrac"
## 
## $DPCoA
## [1] "dpcoa"
## 
## $JSD
## [1] "jsd"
## 
## $vegdist
##  [1] "manhattan"  "euclidean"  "canberra"   "bray"       "kulczynski"
##  [6] "jaccard"    "gower"      "altGower"   "morisita"   "horn"      
## [11] "mountford"  "raup"       "binomial"   "chao"       "cao"       
## 
## $betadiver
##  [1] "w"   "-1"  "c"   "wb"  "r"   "I"   "e"   "t"   "me"  "j"   "sor" "m"  
## [13] "-2"  "co"  "cc"  "g"   "-3"  "l"   "19"  "hk"  "rlb" "sim" "gl"  "z"  
## 
## $dist
## [1] "maximum"   "binary"    "minkowski"
## 
## $designdist
## [1] "ANY"

# Algunos métodos de ORDENACIÓN
Ordenacion <- c("NMDS", "MDS", "PCoA", "DPCoA")
# Algunas DISTANCIAS
Distancia <- c("bray", "unifrac", "dpcoa", "jsd", "jaccard")

#Para poder usar la variable año en la bandera "shape" la convertimos en factor
merged_metagenomes@sam_data@.Data[[6]] <- as.factor(merged_metagenomes@sam_data@.Data[[6]])
for (o in Ordenacion){
  # Prueba 
  # o <- "NMDS"
  for (d in Distancia){
    tryCatch({
      # Prueba 
      # d <- "bray"
      Ord <- ordinate(physeq = merged_metagenomes, method = o, distance = d)
      Beta <- plot_ordination(physeq = merged_metagenomes, ordination = Ord, 
                              color = "SITE_ID", shape = "YEAR") +
        theme_classic() +
        geom_point(size = 5) +
        geom_text(mapping = aes(label = SITE_ID), size = 4, vjust = 2, hjust = 1) + 
        geom_vline(xintercept = 0) +
        geom_hline(yintercept = 0) +
        labs(title = paste("Beta diversity with", o, "-", d, 
                           "for the Calakmul samples ShotGun")) +
               theme(legend.position = "right", 
                     legend.text = element_text(size=12, face="bold"),
                     legend.title = element_text(size = 12, face = "bold"),
                     axis.text = element_text(color = "black", size = 12, face = "bold"),
                     axis.title = element_text(color = "black", size = 12, face = "bold"),
                     plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
               scale_color_brewer(palette = "Set1") + 
               scale_fill_brewer(palette = "Set1")
      print(Beta)
      pdf(paste0("Resultados/BetaDiversidad/PlotBeta", o, "-", d, ".pdf"), 
          width=10, height=7)
      print(Beta)
      dev.off()
      svg(paste0("Resultados/BetaDiversidad/PlotBeta", o, "-", d, ".svg"),
          width=10, height=7)
      print(Beta)
      dev.off()
    }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) 
  }
}

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.02778263  max resid 0.02891375 
## Run 2 stress 0 
## ... Procrustes: rmse 0.3444346  max resid 0.4253598 
## Run 3 stress 0 
## ... Procrustes: rmse 0.06473031  max resid 0.07607818 
## Run 4 stress 0 
## ... Procrustes: rmse 0.02035538  max resid 0.02259239 
## Run 5 stress 0 
## ... Procrustes: rmse 0.06570887  max resid 0.07694945 
## Run 6 stress 0 
## ... Procrustes: rmse 0.1943052  max resid 0.2572746 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1172204  max resid 0.1400096 
## Run 8 stress 0 
## ... Procrustes: rmse 0.05258035  max resid 0.0601127 
## Run 9 stress 0 
## ... Procrustes: rmse 0.3032284  max resid 0.4092741 
## Run 10 stress 0 
## ... Procrustes: rmse 0.07178567  max resid 0.08471621 
## Run 11 stress 0 
## ... Procrustes: rmse 0.3213424  max resid 0.4499814 
## Run 12 stress 0 
## ... Procrustes: rmse 0.3369762  max resid 0.47323 
## Run 13 stress 0 
## ... Procrustes: rmse 0.1224177  max resid 0.1528506 
## Run 14 stress 0 
## ... Procrustes: rmse 0.1830156  max resid 0.2293378 
## Run 15 stress 0 
## ... Procrustes: rmse 0.2391714  max resid 0.3252009 
## Run 16 stress 0 
## ... Procrustes: rmse 0.1224816  max resid 0.1484967 
## Run 17 stress 0 
## ... Procrustes: rmse 0.1958854  max resid 0.2594873 
## Run 18 stress 0 
## ... Procrustes: rmse 0.1352705  max resid 0.1625515 
## Run 19 stress 0 
## ... Procrustes: rmse 0.3467075  max resid 0.4864014 
## Run 20 stress 0 
## ... Procrustes: rmse 0.07094986  max resid 0.0839901 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold

## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.

## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty. 
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.2344238  max resid 0.2937118 
## Run 2 stress 0 
## ... Procrustes: rmse 0.1935433  max resid 0.2591064 
## Run 3 stress 0 
## ... Procrustes: rmse 0.3277101  max resid 0.3423276 
## Run 4 stress 0 
## ... Procrustes: rmse 0.2907152  max resid 0.3136829 
## Run 5 stress 0 
## ... Procrustes: rmse 0.2350237  max resid 0.2776408 
## Run 6 stress 0 
## ... Procrustes: rmse 0.2392581  max resid 0.2671299 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1867016  max resid 0.2618725 
## Run 8 stress 0 
## ... Procrustes: rmse 0.2139299  max resid 0.2524363 
## Run 9 stress 0 
## ... Procrustes: rmse 0.1946998  max resid 0.2310102 
## Run 10 stress 0 
## ... Procrustes: rmse 0.2921048  max resid 0.3224547 
## Run 11 stress 0 
## ... Procrustes: rmse 0.1562071  max resid 0.2136159 
## Run 12 stress 0 
## ... Procrustes: rmse 0.2650681  max resid 0.3206244 
## Run 13 stress 0 
## ... Procrustes: rmse 0.169596  max resid 0.2018941 
## Run 14 stress 0 
## ... Procrustes: rmse 0.2484225  max resid 0.2982738 
## Run 15 stress 0 
## ... Procrustes: rmse 0.04849595  max resid 0.06031097 
## Run 16 stress 0 
## ... Procrustes: rmse 0.1296971  max resid 0.1542298 
## Run 17 stress 0 
## ... Procrustes: rmse 0.2222494  max resid 0.253172 
## Run 18 stress 0 
## ... Procrustes: rmse 0.07582651  max resid 0.09313853 
## Run 19 stress 0 
## ... Procrustes: rmse 0.02437069  max resid 0.03032981 
## Run 20 stress 0 
## ... Procrustes: rmse 0.04372944  max resid 0.05389333 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin

## Warning in metaMDS(ps.dist): stress is (nearly) zero: you may have insufficient
## data

## Warning in metaMDS(ps.dist): Shape variable was not found in the available data
## you provided.No shape mapped.